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Abstract 


This paper provides a solution to the fundamental linear fractional order differential 
equation, namely, c d?x(t)+ax(t)=bu(t) . The impulse response solution is shown to be a series, 

named the E-function, which generalizes the normal exponential function. The E-function provides 
the basis for a <yth order “fractional pole”. Complex plane behavior is elucidated and a simple 
example, the inductor terminated semi-infinite lossy line, is used to demonstrate the theory. 


Introduction 


The problem to be addressed here is the solution of the fractional order differential equation 

D,' x(t) = e d* x (t) = - a x(t) + bu(t) (1) 

where the notation has been defined in Lorenzo and Hartley (1998). Here it will be assumed for 
clarity that the problem starts at t — 0 , which sets c = 0 . It is also assumed that all initial 
conditions, or initialization functions, are zero. Thus we will primarily be concerned with the 
forced response. The initialization response has been addressed in Lorenzo and Hartley (1998). 
Rewriting Equation (1) with these assumptions gives 

0 d?x(t)= -ax(t) + bu(t). (2) 

We will use Laplace transform techniques to simplify the solution of this differential 
equation. In order to do so for this problem, the Laplace transform of the fractional differential is 
required. Using the results given in Oldham and Spanier (1974) or Lorenzo and Hartley (1998), 
and ignoring initialization terms, Equation (2) can be Laplace transformed as 

s“ X(s)=-aX(s) + bU(s) . (3) 

This equation can be rearranged to obtain the system transfer function 


*<£> -gw— L- 

U(s) s q +a 


(4) 
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This is then the transfer function of the fundamental linear fractional order differential equation. 
As such, it contains the fundamental “fractional” pole (to be discussed later) and is the building 
block for more complicated systems, as discussed in the next paragraph. 

Typically, transfer functions are used to study various properties of a particular system. 
Specifically, they can be inverse transformed to obtain the system impulse response, which can 
then be used with the convolution approach to the problem. Generally, if U (s) is given, then the 
product G(s)U(s) can be expanded using partial fractions, and the forced response obtained by 
inverse transforming each term separately. To accomplish these tasks, it is necessary to obtain 
the inverse transform of Equation (4), which is the impulse response of the fundamental fractional 
order system. 

Unfortunately, referring to standard tables of Laplace transforms, such as Erdelyi (1952) 
or Oberhettinger and Badii (1973), the inverse transform of the right side of Equation (4) is only 
known when q — 0.5, q = 1 .0, q = 2.0, or a - 0 . As the intention is to obtain the solution for 
arbitrary q , it is necessary to derive generalized fundamental impulse response for the fractional 
order differential equation, Equation (2). This is done in the next section using Laplace 
transforms. 


The Generalized Impulse Response Function 

Although the Laplace transform tables do not contain terms of the form of Equation (4), 
they do contain the transform pair 



q >0. 


(5) 


Thus, if we can expand the right side of Equation (4) in descending powers of 5, we can then 
inverse transform the series term by term and obtain the generalized impulse response. It should 
be noted that throughout this paper, it is assumed that q > 0 . 

As the constant b in Equation (4) is a constant multiplier, it can be assumed, with no loss 
of generality, to be unity. Then expanding the right side of Equation (4) about s - using long 
division, gives 


G(s) = 


s { + a 


a a 

2a 3 a 

S H S 1 


= 1 y (-«r 

Zj nq 

* «=0 S 


( 6 ) 


This series can now be inverse transformed term by term using Equation (5). The result is 
r 1 a a 2 ] t’ 1 '' 


g(t) = V 


a a 

1 — -f 

s* s 2<l s* q 


ar«-' | aj_ 

H q) T{2q) + r(3q) 


(7) 
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The right side can now be collected into a summation and used as the definition of the 
generalized impulse response function 


oo / \ n f m, 

*<') = ' , ‘ l r«F, [-«.»]. 9>°- 

“ o r(nq + q) 

We also have the important Laplace transform identity 


(8) 


l{f, M}= -r— . 9>0. 

' c- ' n 


(9) 


s' —a 


Here we have defined the notation for this function to be as it is closely related to the 

Mittag-Leffler function E q [at q ], function (Mittag-Leffler, 1903a; Mittag-Leffler, 1903b; 
Mittag-Leffler, 1905). The Mittag-Leffler function is defined as 


*,M« X: 


to nnq + D 

(Erdelyi, 1954). Letting x =-at q , this becomes 


q > 0, 


(10a) 






« f nij 


£Zr(nq + l) 


q> 0 


(10b) 


(Bagley and Calico, 1991), which is similar to, but not quite the same as Equation (8). The 
Laplace transform of this Mittag-Leffler function can also be obtained via term by term transform 
of series (10b), that is 


tfe[- a ,''D=4-L- 


at < 


a 2 t 2q 


1 2 

a a 
+ 


ro) ru+<?) ra+2^) 


5 s q+] s 29+1 


( 11 ) 


or, equivalently 

rk[-»'"]}=- 


. 1 a 

1 H r— + •’• 


s* s lq 


1 c ’ 

/ V' 

— CL y 

S n=0 ' 

\ s / 


( 12 ) 


It should now be recognized that the summation in this expression is similar to Equation (6). 
Using that. Equation (12) can be written as 




(13) 


s +a 
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or equivalently 


Z,{£,[-«r’]}= - [s' £{£,[-«,»]}] (14) 

Thus, the general result can be written 

ifokj'-r- • q>0 . (15) 

s-a 

Also, notice from Equation (14), that the E -function and the F -function can be related as follows, 

. (16) 

This section has shown that the F -function is the impulse response of the fundamental 
linear fractional differential equation. Plots of the F -function and the E -function for various 
values of q are given in Figures I and 2, respectively. 



Time 


Figure 1. The F q [— l,f]-function vs. time as q varies from 0.25 to 2.0 in 0.25 increments. 
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It should be noted that other authors have obtained a solution to Equation (2), but they appear to 
be much less direct. Bagley and Calico (1991) obtain a solution in terms of Mittag-Leffler 
functions. Miller and Ross (1993) obtain a solution in terms of the fractional derivative of the 
exponential function. They use the function 

E,(v,a)= 0 d;"e“. (17) 

whose Laplace transform is 

L{E>, 0 )}=— . ( 18 ) 

s-a 

Also, Glockle and Nonnenmacher (1991) obtain a solution in terms of the even more complicated 
Fox Functions. Clearly, all of these functions are useful for this problem (Eqn. (4)), but the 
F -function presented here appears to most properly generalize the exponential function for use 
with fractional differential equations. Finally, it should be noted that the F -function is also 
mentioned by other authors as well. Oldham and Spanier (1974, page 122) mention it in passing 
in a footnote discussing eigenfunctions. We have recently discovered that Robotnov (1980) and 
(1969) studied the F -function extensively with respect to hereditary integrals (he calls it the 
Cyrillic backwards E, or “eh”-function) . Our assumption is that the fractional calculus 
community has not discovered this work as it has been “hidden” there. In the next section we 
will consider various properties of this function. 



Figure 2. The Mittag-Leffler function, E q [-t q ] vs. time as q varies 
from 0.25 to 2.0 in increments of 0.25. 
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Important Properties of the F-function 


In this section, several properties of the F -function are derived. This is done specifically, 
so that it can be shown that the F -function solves the fractional differential equation, Equation (2), 
by direct substitution. It is also shown that the F -function satisfies what Oldham and Spanier 
(1974) refer to as the “eigenfunction” property. This essentially means that the q' h -derivative of 
the function F q [a, t 11 ], returns the same function F q [a, t q ] for t > 0, (see Equation (27)). Several 
intermediate results are necessary to show these properties, and they are now derived. 

First of all, we will consider the step response of the system given in Equation (2). This 
can be obtained via Laplace transforms by transforming the input function u(t), which is chosen 
to be a unit step function. Its Laplace transform is Y , which must then be multiplied by the 
transfer function, via Equation (4), to give the transformed step response as 


*(*) = - 

s 


1 

s'' +a 


(19) 


This equation can be manipulated to give 


*(*) = 


1 la 

a 

1 / a 

1 _ * 

1 

: a 

l 

' 

s 

s q +a 

s 

s“+a_ 

s 


( 20 ) 


This equation can now be inverse transformed using Equation (14) for the second term on the 
right. The result is the step response of the system. 


a [s^+a) 


( 21 ) 


where H(t) is the Heaviside unit step function, and which also gives another Laplace transform 
identity. This step response is given in Figure 3 for several values of q and a = 1. It is also 
interesting to notice that taking the integer derivative ( 0 r/,') of both sides of Equation (21) 

necessarily gives the F -function on the left (the derivative of the step response is the impulse 
response), and a new identity on the right; 

Fj-a,t]=±^-(H(t)-E i/ [-at‘']). (22) 

a at 


Now referring back to Equation (14), and multiplying the Laplace transform there by s Ogives 

s~ lf l{e \~at q }=*■' L{F\-a,t]}=-J- — t , which is the Laplace transform of Equation (19). 

' ,y(s’+aj 
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Inverse transforming this using Equation (21) shows that the step response of the system of 
Equation (2) is also equal to the q th -integral of the Mittag-Leffler function; that is 


r'j-J— .I = [-«*]= 0 rf;’ £,[-«''] . (23) 

[ .v [s { + a ) j a 

Some other interesting identities can be obtained by taking the q th -derivative of the 
F -function. Taking the uninitialized derivative ( 0 d q ) in the Laplace domain by multiplying 
by s q gives 


r 1 





(24) 


It should now be noticed that this is also the integer derivative of the Mittag-Leffler function, 


L~' 





(25) 



Time 


Figure 3. The step response of the system of Equation (2) as q varies 
from 0.25 to 2.0 in 0.25 increments. 
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This equation can also be rewritten as 


\ = 8(t)-aF 9 [-a,t] , (26) 

[s" + a] { s q + a J q 

where the delta function is recognized as the unit impulse. Now comparing Equations (24) and 
(26), it can be seen that 

of/," F (l [-a,t]= 8(t) - aF q [-a,t] . (27) 

This equation demonstrates the eigenfunction property of returning the same function upon 
q h -order differentiation. This is a generalization of the exponential function in integer order 
calculus. 

It is now easy to show that the F -function is indeed the impulse response of the system 
of Equation (2). Referring back to Equation (2), inserting u(t) = S(t), and setting b- 1 , yields 

o d* x(t) = -ax{t) + <5(0« (28) 

For the F -function to be the impulse response of the system, it must be the solution to 
Equation (2), that is x(t) = F Inserting this into Equation (28) gives 

0 <i;' F q [-a,t]= -aF ll [-a,t]+ 8{t) . (29) 

This equation has been obtained by direct substitution into the differential equation. Referring 
back to Equation (27), however, shows that the q th -derivative of the F -function on the left is in 
fact equal to the right side of Equation (29). Thus it is shown by direct substitution that the 
F -function is indeed the impulse response of the system of Equation (2). 


Behavior of the F-function as the Parameter a Varies 

In 1903, Mittag-Leffler (1903a) introduced his new function E q \ax\ q > 0 . He considered the 

parameter a to be a complex number, a . As he studied this function (Mittag-Leffler, 

1903b; Mittag-Leffler, 1905) it became apparent that this function was either stable (decays to 
zero) or unstable (goes to infinity) as x increases, depending upon how he chose the parameters 
a and q . The result was that the function remained bounded for increasing x if 

• < 30 > 

This section will demonstrate that the F -function shares this property. Furthermore, this result 
carries over directly to the Laplace s-domain, which provides a fairly straightforward approach 
for proving the result of Mittag-Leffler in Equation (30). 
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The discussion follows most easily if we consider the Laplace transform representation of 
the F -function from above, 

^[- 0 , 1 ]}=-^—, q> 0. (3D 

s' + a 

Generally, to understand the dynamics of any particular system, we often consider the nature of 
the s-domain singularities. We will define s = re 1 6 in what follows. The particular function of 
Equation (31), however, does not have any poles on the primary Riemann sheet of the s-plane 
(|e|<7T), as it is impossible to force the denominator of the right side of Equation (31) to zero 

anywhere in the s-plane. Notice, however, that it is possible to force the denominator to zero if 
secondary Riemann sheets are considered. For example, the denominator of the Laplace 
transform 



(32) 


does not go to zero anywhere on the primary sheet of the s-plane, (|0|<7T ). It does go to zero on 
the secondary sheet, however. With s=e ±j2n , the denominator is indeed zero. Thus this 
Laplace transform has a pole at s=e ±j2 * , which is at s - 1 + jO on the second Riemann sheet. 


This is shown in Figure 4, where 


r* + l 


is plotted as a function of Real(s) and Imaginary(s). 


Normally, to get to a secondary Riemann sheet, it is necessary to go through a branch cut 
on the primary Riemann sheet. This is accomplished by increasing the angle in the s-plane. 
Referring to Figure 4, increasing the angle until 0=+7t , gets us to the branch cut on the s-plane. 
This can also be accomplished by decreasing the angle until 6=-n , which also gets us to the 
branch cut. Thus the branch cut lies at s = re ±> * , for all positive r . Further increasing the angle 
eventually gets us to 6 =± 2lZ . For the Laplace transform of Equation (32), the behavior of the 
transform is completely described with the two Riemann sheets. Returning to the primary 
Riemann sheet on the s-plane, the branch cut begins at 5 = 0 , the s-plane origin, and extends out 
the negative real axis to infinity. The ends of the branch cut are called branch points, which are 
then at the origin and at minus infinity in the s-plane, for the example. These branch points can 
also be considered to be singularities on the primary sheet of the s-plane as well, but the Laplace 
domain function does not go to infinity there. Whereas inverse Laplace transforms are usually 
obtained by integrating around these branches and branch points on the primary sheet, this 
thinking effectively ignores the secondary sheets, where the singularities (poles) are located. 
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Figure 4. Both sheets of the Laplace transform of the F -function in the s-plane. 


As it is difficult to visualize multiple Riemann sheets, it is useful to perform a conformal 
transformation, following LePage (1961), into a new plane. Here we will let 

W5 5*. (33) 

The transform in Equation (31) then becomes 






s H +a 


w+a 


(34) 


With this transformation, we will study the w-plane poles. Once we understand the time domain 
responses that correspond to the w-plane pole locations, we will be able to clearly understand the 
implications of this new complex plane. 

To accomplish this, it is necessary to map the s-plane, along with the time-domain 
function properties associated with each point, into the new complex w-plane. To simplify 
discussion we will limit the order of the fractional operator to 0 < q < 1 . Let 

w= pe J * = (X+jp . (35) 
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Then referring to Equation (33) 


w = s“ = (re' 9 ) 1 = r“e iq6 . (36) 

With this equation, it is possible to map either lines of constant radius, or lines of constant angle 
from the s-plane into the w-plane. Of particular interest, is the image of the line of s-plane 

stability (the imaginary axis), that is s = re ±j ' 2 . The image of this line in the w-plane is 


w = r'‘e +i *^ 


(37) 


which is the pair of lines at (p=—~ . Thus, the right half of the s-plane maps into a wedge in 
the w-plane of angle less than ±90q degrees, that is, the right half s-plane maps into 



(38) 


For example, with q = / 2 , the right half of the s-plane maps into the wedge bounded by 
4>< q *A , see Figure 5. 

It is also important to consider the mapping of the negative real s-plane axis, 5 = re ,K . 
The image is 


w = r q e ±)q * . (39) 

Thus the entire primary sheet of the s-plane maps into a w-plane wedge of angle less than 
±180<7 degrees. For example, if q = / 2 , then the negative real s-plane axis maps into the 

w-plane lines at ± 90 degrees, see Figure 5. 

Continuing with the q = / 2 example, and referring to Figure 5, it should now be clear that 
the right half of the w-plane corresponds to the primary sheet of the Laplace s-plane. All of the 
time responses we are familiar with from integer order systems have poles that are in the right 
half of the w-plane. The left half of the w-plane however, corresponds to the secondary Riemann 
sheet of the s-plane. A pole at w = -l + jO lies at s = +\ + jO, on the secondary Riemann 
sheet of the s-plane. This point in the s-plane is really not in the right half s-plane, corresponding 
to instability, but rather is “underneath” the primary s-plane Riemann sheet, or even more 
intuitively satisfying, this point is “inside” the negative real s-plane axis. Lying inside the 
negative real s-plane axis is a better image, as the easiest way to get to this pole is by increasing 
the angle of an s-plane point until d=±7t , at which time you are on the negative real s-axis. 
Increasing the angle any further takes you “inside” the negative real s-axis onto the secondary 
Riemann sheet, and consequently farther away from the right half s-plane. As the corresponding 
time responses must then be even more than over-damped, we will call any time response whose 
pole is on a secondary Riemann sheet, “hyperdamped.” It should now be easy for the reader to 
extend this analysis to other values of q . 
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Figure 5. w-plane for q = Y 2 . 


To summarize the above, the shape of the F -function time response, F [— aj], depends 
upon both q , and the parameter - a , which is the pole of Equation (34). This is shown in 
Figure 6. For a fixed value of q , the angle (f> of the parameter -a , as measured from the positive 
real w-axis, determines the type of response to expect. For small angles, the time 

response will be unstable and oscillatory, corresponding to poles in the right half s-plane. For 
larger angles, the time response will be stable and oscillatory, corresponding to 

poles in the left half s-plane. For even larger angles, |0|>g;r, the time response will be 
hyperdamped, corresponding to poles on secondary Riemann sheets. 

It is now possible to do fractional system analysis and design directly in the w-plane. To 
do this, it is necessary to first choose the greatest common fraction (q) of a particular system 
(clearly non-rationally related powers are a problem and will be considered in a future paper, 
although a close approximation of the irrational number will be sufficient for practical 

application). Once this is done, all powers of s q are replaced by powers of w . Then the standard 
pole-zero analysis procedures can be done with the w-variable, being careful to recognize the 
different areas of the particular w-plane. This analysis includes root finding, partial fractions 
(note that complex conjugate w-plane poles still occur in pairs), root locus, compensation, etc. 
We have thus now completely characterized all possible behaviors for fractional order systems in 
a new complex w-plane; that is, given a set of w-plane poles, the corresponding time domain 
functions are known both quantitatively and qualitatively. Although most of the discussion has 
actually been for 0 < q < J4 , it is somewhat applicable to larger values of q with the appropriate 
modifications for many-to-many mappings. 
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Example 


In this section, a simple example is presented to demonstrate how to use the theory 
presented in this paper to obtain the solution to a physical fractional order system. The system 
considered is the inductor terminated lossy line studied by Heaviside (1922) and Bush (1929) and 
shown in Figure 7. 


sL 


r-Yl 


&) V M 


+ wWVVtWVi — 

v a (s) 4 = 4 = 4 = ■=> 


Figure 7. Inductor terminated semi-infinite lossy line example. 


The input to the system is the voltage at the left, and the output will be chosen to be the voltage at 
the terminal of the lossy line. Using impedances, with L=l, gives the transfer relationship to be 


V,(s) 


= G(s) = 


77 


s + 


77 


s*+ 1 


(40a) 


It should be noted that this problem can be written in the time domain as 


„d% vjt) + v„(r) = v,(/) 


(40b) 


where it is assumed that the initializations are zero. 


This problem can be solved in several ways depending upon the specific input and also 
depending upon the base value of q that is chosen. Clearly, with q = Y » the impulse response of 
the system is given by Equation (9) as 


v 0 {t) = L 


,X 


+ 1 


= Fy\-U] 


(41) 


The shape of this function can be seen in Figure 1 . The step response of this system can also be 
found from Equation (21) as 


v 0 (t) = L-' 



= «</)-£*[-»*]. 


(42) 


The shape of this function can be seen in Figure 3. 
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It is now instructive to solve this problem by assuming that the basis q = x /i instead of 
q = y 2 . Clearly the answer must be the same; however, this approach will demonstrate the use of 
the w-plane on a system with a known response, as well as some other interesting properties of 
the F -function. In the original transfer function in Equation (40), let s = w^. The transfer 
function then becomes 


VpM 

V,(w) 


= G(w) = 



(43) 


This transfer function has w-plane poles at w — —l, w — e J// ' , w — e y/ ' . Referring back to the 
w-plane in Figure 5, all the w-plane poles are to the left of the instability wedge at <p=±Y A . The 

two poles in the right half w-plane correspond to s-plane poles at s = e J , and thus an 

oscillatory response is expected. The third pole at w—~ 1, is in the hyperdamped region and 
should indicate a rapidly decaying time response added to the oscillatory response from above. 
To obtain the impulse response, the w-plane transfer function must be expanded in partial 
fractions using the base value of q , 


G(w) = 


0.3333 
w + 1 


0.1667 + j0.2887 
w - 0.5000 - 79.8660 


0.1667 - 70.2887 
w - 0.5000 + 70.8660 ' 


(44) 


The corresponding time response can be obtained by inverse transforming term by term to give 

v 0 (t) = 0.3333 Fy\-\,t]- (0.1667 + j0.2887)fJo. 5000+ 7 0.8660,f] 

- (0.1667 - 70.2887)/^ [0.5000- 70.8660,r] ' 

which is equivalent to the solution given in Equation (41). This results also demonstrates that 
F -functions of different indices can be directly related to one another. 


Summary 

The fundamental linear fractional order differential equation has been considered and its 
impulse response has been obtained as an F -function. Several properties of this function have 
been presented and discussed. In particular, the Laplace transform properties of the F -function 
have been discussed using multiple Riemann sheets and a conformal mapping into a more readily 
useful complex w-plane. 

It is felt that this generalization of the exponential function, the F -function, is the most 
easily understood and most readily implemented of the several other generalizations presented in 
the literature. 
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